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Many-body systems with both coherent dynamics and dissipation constitute a rich class of models 
which are nevertheless much less explored than their dissipationless counterparts. The advent of 
numerous experimental platforms that simulate such dynamics poses an immediate challenge to sys¬ 
tematically understand and classify these models. In particular, nontrivial many-body states emerge 
as steady states under non-equilibrium dynamics. While these states and their phase transitions 
have been studied extensively with mean field theory, the validity of the mean field approximation 
has not been systematically investigated. In this paper, we employ a field-theoretic approach based 
on the Keldysh formalism to study nonequilibrium phases and phase transitions in a variety of mod¬ 
els. In all cases, a complete description via the Keldysh formalism indicates a partial or complete 
failure of the mean field analysis. Furthermore, we find that an effective temperature emerges as 
a result of dissipation, and the universal behavior including the dynamics near the steady state is 
generically described by a thermodynamic universality class. 

PACS numbers: 64.60.Ht, 64.70.qj 


I. INTRODUCTION 

Condensed matter systems typically relax to their equi¬ 
librium state on very short time scales. A success¬ 
ful paradigm of statistical physics developed over many 
years explains various aspects of equilibrium or near¬ 
equilibrium phenomena in many-body systems. Specifi¬ 
cally, quantum phase transitions, which emerge at zero 
temperature, have been the subject of intense research 
over the past several decades [1]. 

In contrast, experiments with ultracold matter have 
opened new avenues to probe far-from-equilibrium many- 
body systems in the presence of both coherent dy¬ 
namics and controlled dissipation, the so-called driven- 
dissipative systems. While the subject is not new, 
the vast range of experimental platforms have brought 
driven-dissipative many-body systems into light again 
and made it necessary to undertake a more thorough 
investigation. Experiments range from polariton con¬ 
densates in the context of semiconductor quantum wells 
in optical cavities [2-6] and arrays of microcavities [7- 
10] to trapped ions [11, 12] and optomechanical se¬ 
tups [13-15]. Furthermore, experiments on strongly in¬ 
teracting Rydberg polaritons are already probing non¬ 
equilibrium many-body physics [16-19], while cavity- 
quantum-electrodynamics experiments can potentially 
explore dynamics under exotic many-body Hamiltonians 
with glassy ground states [20-22]. The interplay be¬ 
tween dissipation, which is generically present in these 
systems, and coherent dynamics leads to new, and in¬ 
herently nonequilibrium, phases. The fundamental ques¬ 
tion is then how we should understand and classify these 
phases. Do dissipative-driven systems give rise to new 
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phases of matter? What are the universal properties of 
the associated phase transitions? 

The sophisticated toolbox of equilibrium physics is 
not immediately applicable to nonequilibrium problems. 
The equivalents of equilibrium concepts such as tem¬ 
perature, free energy, and partition function either have 
no obvious counterpart out of equilibrium, or often be¬ 
come intractably complicated. In particular, the main 
goal of condensed matter physics is to find properties 
of the ground state or, at finite temperature, the ther¬ 
mal state of the system, while, out of equilibrium, non¬ 
trivial many-body states emerge as steady states under 
non-equilibrium dynamics. In the absence of a powerful 
systematic approach, approximations such as mean field 
theory have been widely used [23-31], but are often at 
odds with other analytical and numerical studies [32, 33], 
sometimes even in the limit of infinite dimensions [34] 
where the equilibrium mean field is known to be exact; 
these analytical and numerical studies are based on vari¬ 
ational methods [32, 33] and approximate rate equations 
[34]. In general, numerical techniques are either limited 
to one dimension, e.g. t-DMRG [35], or to infinite dimen¬ 
sions such as nonequilibrium dynamical mean field theory 
[36, 37], or cannot be applied to non-equilibrium systems, 
e.g. path-integral Monte Garlo [38]. It is worth mention¬ 
ing that there exist exact solutions, due to integrability, 
for a very specific class of nonequilibrium models where 
the system is driven only at the boundaries [39, 40]. 

In principle, the Keldysh-Schwinger functional integral 
formalism provides a general approach to nonequilibrium 
physics. A notable example is the universal behavior of 
early evolution of an initial state prepared far from equi¬ 
librium [41], see also [42] for a review. This approach has 
been applied to a number of driven-dissipative systems 
such as lossy polariton condensates [43-45] and driven 
atomic ensembles interacting with a cavity mode [46]. 
Indeed, a systematic application of the Keldysh formal- 
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Model 

Mean Field 

Field Theory 

A 

H = ^ Ei erf 

Li = \/T a~ 

Continuous transition from 
an ordered phase ((erf) =const) to 
a disordered phase ((erf) = 0) 

First order transition 

at sufficiently strong dissipation 

B 

H = - JEfe-> H E, + A E, 

Li = vT o~ 

Continuous transition from 

a phase with one stable solution to 
a bistable regime 

No bistability 

C 

Li = VT G~ 

Continuous transition from 
an XY phase (staggered {of)) to 
a disordered phase {{of) = 0) 

No XY phase in 

d = 2 dimensions 

D 

H = -JY.{ij)^t<^J TAEiCrf 

L\ = y/fa- Lf = yr;<T+ 

Lij = (1 -I- aj) 

No phase transition 
or spontaneous symmetry breaking 
{{Oi) = {oD = 0) 

A continuous transition to a 
spontaneously symmetry broken 
phase {{of) = const) 


Table I: Summary of the results. We consider four driven-dissipative models with spin-l/2s on a d-dimensional cubic lattice. 
In each case, the Hamiltonian (H) and the dissipation via the Lindblad operators (L) are defined. The mean-field theory 
prediction is described and contrasted with the Keldysh field-theoretic treatment. Nearest-neighbor Ising, isotropic XY, and 
anisotropic XY interactions are considered in different examples. In models (A), (B), and (C), we take d = 2 or 3, while, in 
model (D), we take d = 3. In the former three models, the dissipation is an independent spontaneous emission on each site, 
whereas, in the latter model, pump and correlated dissipation on nearby sites are assumed as well. The interplay of unitary 
and dissipative dynamics creates rich phase diagrams. Our analysis indicates that the mean-held predictions are incomplete or 
even wrong for all four models, and should be supplemented with held theory via the Keldysh formalism. 


ism to the wide variety of driven-dissipative problems is 
far from complete. Surprisingly, even the simplest many- 
body driven-dissipative spin models have not been fully 
tackled with the Keldysh formalism; these models should 
not be confused with spin-boson models where a two-level 
system is coupled to a dissipative environment in ther¬ 
mal equilibrium [47]. To be more precise, in contrast to 
spin-boson models where spins are strongly coupled to a 
thermal environment, driven-dissipative systems studied 
in this paper deal with situations where coupling to the 
environment is weak, but the system is driven out of ther¬ 
mal equilibrium by an external drive (for example, by a 
laser beam), or it may be in contact with two different 
baths at different temperatures, or may be coupled to a 
non-thermal bath. 

In this paper, we consider a variety of nonequilibrium 
models, mostly spin models on a d-dimensional cubic lat¬ 
tice, and study their, inherently nonequilibrium, steady 
states. The fact that we are interested in the near-critical 
(criticality identified by a diverging time scale) long¬ 


distance behavior of the models allows us to map the spin 
models to continuum field theories, which we study via 
the Keldysh formalism in great detail. For each model, 
we compare and contrast the field-theoretic Keldysh ap¬ 
proach to the mean field theory which is shown to miss 
some or most of the features of the full field-theory treat¬ 
ment. A unified and systematic Keldysh approach is ap¬ 
plied to most of these models (Models A-C in Table I): 
Close to the mean-field phase transition, critical and mas¬ 
sive components of the field are identified, and the latter 
is integrated out to find an effective action for the critical 
field. While such a procedure for obtaining an effective 
action is standard, the Keldysh formalism involves two 
different, classical and quantum, fields which require spe¬ 
cial care. 

We first provide a brief description of the tools used 
throughout this paper in Sec. II, and then study four spe¬ 
cific driven-dissipative models in Sec. III. We generically 
find that a driven-dissipative model behaves thermally, 
that is, an effective temperature emerges, and the phase 
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transition between different phases is described by well- 
known thermodynamic paradigms and their universality 
classes. The emergence of an effective equilibrium and a 
conventional thermal phase transition has been identified 
in many contexts [21, 25, 43, 44, 46, 48-53]. Further¬ 
more, effects beyond mean field have been studied in the 
context of driven-dissipative condensates [43-45] . Never¬ 
theless, more generally, the determination of what precise 
thermodynamic phase corresponds to the nonequilibrium 
steady state is often nontrivial and depends sensitively 
on the type of dissipation and its competition with the 
coherent dynamics. Systematic derivation of this cor¬ 
respondence for each of the models under consideration 
constitutes the main result of the present manuscript. 
What is universal in all these models is the emergence 
of thermal phase transitions, but they show generic, al¬ 
though model-dependent, effects of fluctuations in driv¬ 
ing phase transitions to a different order (model A), re¬ 
moving mean-field artifacts (models B, C), melting order 
(model C), or inducing symmetry breaking (model D). 

For the benefit of the reader, we summarize our results 
in Table I, which introduces the four driven-dissipative 
models under consideration and highlights the important 
differences between the results obtained via mean field 
theory and via the Keldysh field-theoretic formalism. 


II. GENERAL FRAMEWORK 


In this section, we closely follow Refs. [23, 43, 44, 46] to 
introduce the general framework in which we define and 
study driven-dissipative systems; see also the pioneering 
works in Refs. [54-56]. To properly treat dissipation, we 
first introduce a second-quantized master equation under 
which the density matrix evolves with both unitary and 
dissipative dynamics. We briefly sketch how the master 
equation is mapped to the Keldysh path integral. We 
then consider a generic form of the Keldysh action, and 
present a simple scaling argument that constitutes the 
first step of the RG procedure. Finally, we show that, 
under certain conditions, the Keldysh path integral maps 
to a classical Langevin equation, and the resulting steady 
state can be expressed as a classical partition function. 
While the tools and methods discussed in this section are 
known and have been utilized in the literature, we find it 
useful to present them in an introductory section rather 
than an appendix as we will use them frequently, and 
mostly in the context of models to which they have not 
been directly applied. For a detailed introduction to the 
Keldsyh formalism, we refer the reader to Refs. [57, 58]. 
The application of the Keldysh formalism to a number of 
driven-dissipative systems can be found in Refs. [43, 44, 
46] , which this work has greatly benefitted from. 

Quantum master equation .—To describe an open sys¬ 
tem, one should include both the coherent and dissipative 
processes in a master equation for the density matrix as 


(in units where = 1) 

dtp = -i[H, p] + ^Yl “ L^LaP - pL^La) ■ 

a 

( 1 ) 

The first term on the right-hand side gives the usual co¬ 
herent evolution via the Hamiltonian H. The dissipation 
is subsumed in the second term (sometimes referred to 
as the linear Liouville operator £[•] acting on p) char¬ 
acterized by the so-called Lindblad operators Lq, which 
describe an incoherent process a. The master equation 
relies on the Born-Markov approximation, which assumes 
that the reservoir has a short correlation time and is large 
enough to be unaffected by the coupling to the system 
[59]. In the so-called stochastic-wavefunction interpreta¬ 
tion, the operator acts as an occasional, discontin¬ 
uous, jump from one state to another as a result of an 
incoherent process [60] (more precisely, the first term in 
parentheses in Eq. (I) describes such a jump, while the 
last two terms take dissipation into account between the 
jumps). Typical Lindblad operators are local operators 
causing transitions between levels or decay of coherences, 
i.e. dephasing. Various examples of Hamiltonian dynam¬ 
ics and Lindblad operators are discussed in Sec. HI and 
summarized in Table 1. 

Keldysh functional integral .—The Keldysh formalism 
provides a general framework to study nonequilibrium 
problems with functional integral techniques. Within 
this approach, the action is defined on a closed time con¬ 
tour with a forward and a backward branch. In perform¬ 
ing the functional integration, one should sum over all 
configurations with independent values of the underlying 
fields on the two branches. To be concrete, consider {a} 
as a shorthand for all the fields in a particular model. 
The Keldysh path integral can be formally expressed as 

where cr = ± denotes the forward and backward 
branches, respectively, while (t, x) are the time and spa¬ 
tial coordinates. Evaluating the functional integral on a 
closed time contour means that the values of a± should 
match at t = oo (at t = — oo, the system is described by 
an initial state whose precise form is unimportant for the 
steady state of the system at long times [57], although 
the early evolution of the system depends sensitively, and 
perhaps even universally, on the initial state [41]). All the 
information in the, possibly time-dependent, state and 
specifically all correlation functions can be computed by 
inserting fields in the functional integration; for a general 
operator O, one has 


Tr 


p{t)0 


(o+W), 


( 3 ) 


where the expectation value (•) is computed with respect 
to the functional integral in Eq. (2). The subscript -|- on 
O on the right-hand side of the above equation implies 
that all the fields inside O are evaluated on the forward 
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branch; one can equally well arrange to have the observ¬ 
able on the backward branch, but it is, in fact, often 
more convenient to work in a different basis, which we 
shall define shortly [57, 58]. 

Conveniently, the master equation (1) can be directly 
mapped to a Keldysh action comprising both the coher¬ 
ent {H) and dissipative (D) terms as 


Sk = Sh + Sd ■ (4) 

The coherent part of the action can be written in the 
coherent-state representation (after normal ordering, and 
assuming that a corresponds to a bosonic operator) of the 
path integral as 


Sh= ^ 

O' —+ ,— 

The relative sign of the forward and backward branches 
has its origin in the minus sign in the commutator [77, p\. 
The dissipative dynamics yields the Keldysh term 

sd = -^Y.I \ , 

( 6 ) 

with the Lindblad terms La and L* given in terms of 
position- and time-dependent fields (see Ref. [44] for 
more details). In general, operators acting on the density 
matrix p from the left (right) give rise to a corresponding 
term on the cr = -I- (cr = —) contour [44, 46]. 

It is often more convenient to work in the Keldysh basis 
defined as 


alidttta I - 77[a<^,a*] 


( 5 ) 


_ a+ -b a_ _ a+ - a_ 

aez- v/^ , a,- v/_ , (7) 

where a^i/q are the so-called classical and quantum fields; 
typically, in an ordered phase {ad) = const, while (a^) = 
0, that is, the quantum field is purely fluctuating. 

Scaling dimensions .—We often encounter the Keldysh 
action of the form ^ a*{idt -b — i^)aci + iTjag]^ at 
the quadratic level, where T and T' generally designate 
dissipation rates in the model. To find the scaling dimen¬ 
sion of the fields as a first step of the RG procedure, we 
perform a simple dimensional analysis; setting the scal¬ 
ing dimension of spatial and time derivatives as [V] = 1 
and [9t] = z, we have z = 2 a.t the quadratic order. We 
also choose dim[r] = 0. The scaling dimensions of the 
classical and quantum fields are then [43, 44] 

r 1 d — 2 d -b 2 

[aci\ = , [og] = , (8) 

with d the number of spatial dimensions. Notice that, 
with our choice of renormalization, the term ^ 
only marginally relevant; any additional powers of Od/q 
or higher-order derivatives in the integrand make it ir¬ 
relevant [102]. Therefore, the quantum held Ug appears 
at most quadratically, and without any derivatives at the 


quadratic order, in the relevant part of the action. We 
stress that our analysis based on power counting is valid 
at or near criticality (criticality dehned as T' = 0). 

Langevin equation .—In the models presented in this 
paper, we often Hnd a Keldysh path integral of the form 


D[ad/q{t,x.)] e 


[tlci dg] 


with 



^ad 


Sfjaci] 

Ki J 


-b c.c. -b iTjog 


(9) 


where the quantum held appears at most at the quadratic 
level, and / is assumed to be a complex-valued functional 
of Qd ; / may generically include gradients of the held as 
well as interaction terms. By the virtue of a Hubbard- 
Stratonovich transformation, one can cast the quadratic 
term in the quantum held in terms of a ‘noise’ held ^(t, x). 
The huctuations of the classical held in the above path 
integral can be then cast exactly in terms of a classical 
Langevin equation as [57] 


idtad{t,x.) = +C(i,x), 

where noise correlations satisfy 


(^(^,x)^*(t',x')) = r(5(x - x)d{t - t'). 


( 10 ) 


( 11 ) 


All correlations of the classical held can be computed 
either via the Keldysh path integral, Eq. (9), or equiva¬ 
lently via the Langevin equation (10). The latter, how¬ 
ever, has the obvious advantage of mapping to a classical 
stochastic equation. We are often interested in hnding 
the steady state of Eq. (10), which, via a series of steps 
outlined above, is indeed equivalent (for the purposes of 
probing the long-wavelength physics near criticality) to 
the steady state of the quantum master equation (1). The 
calculation of the steady state of Eq. (10) is not trivial in 
general. However, for the special case—encountered in 
Models A-C in Table I — where the function / is purely 
imaginary [103], f[ad] = ifi[ad], the steady state is given 
by the probability distribution 

P[ad] ~ exp , (12) 

which takes the form of a thermodynamic partition func¬ 
tion with the effective temperature = r/2, and // 
as the effective classical Hamiltonian of the system. The 
value of the effective temperature is chosen such that 
the fluctuation-dissipation condition, relating the fluc¬ 
tuations of the noise to the dissipation in the effective 
model, holds. Computing correlation functions weighted 
by Eq. (12) will produce all the relevant information in 
the quantum steady state. We also remark that Eq. (12) 
can be derived from the Fokker-Planck equation that 
casts the Langevin equation (10) in the form of an equa¬ 
tion for the probability distribution [61]. 

We stress that driven-dissipative systems may not al¬ 
ways be described by a thermal-like probability distribu¬ 
tion. In particular, this description may fail if the system 
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is not invariant under a symmetry transformation that 
characterizes the equilibrium condition [44, 53]. 

Although the procedure and the steps outlined in this 
section are rather standard, the key challenge lies in 
bringing a many-body model of interest, e.g. a spin 
model, into the form given by Eq. (9) or Eq. (10); in 
almost all the models studied in this paper (Models A- 
C in Table I; Model D requires a special treatment), we 
find that such a transformation is nontrivial. Neverthe¬ 
less, we show how to achieve this goal via a systematic 
approach. 


III. DRIVEN-DISSIPATIVE MODELS 

Here, we consider a number of driven-dissipative mod¬ 
els, mostly spin systems on a cubic lattice, each defined 
by a particular Hamiltonian and a particular form of dis¬ 
sipation. Several models considered here are novel and 
have not been studied in the literature even at the mean- 
field level. In each subsection, we compare the results of 
the mean-field analysis to those of the full field-theoretic 
treatment using the Keldysh formalism. The mean-field 
analysis is shown to partially or completely fail in all 
the models, while the Keldysh approach provides a field- 
theoretic formalism to go beyond mean field, best suited 
to study the vicinity of phase transitions. 

The driven-dissipative spin models described below can 
be implemented using a variety of experimental systems. 
The Hamiltonians can be implemented using ions cou¬ 
pled via motional modes [62-66], atoms in optical lat¬ 
tices coupled via superexchange [67-69] , atoms in optical 
cavities or along waveguides coupled via optical modes 
[21, 70, 71], or polar molecules [72-77], Rydberg atoms 
[78], magnetic atoms [79], magnetic defects in solids [80], 
and Rydberg polaritons [81] coupled via dipole-dipole or 
van der Waals interactions. The Lindblad operators ei¬ 
ther arise in these models naturally via processes like 
spontaneous emission or can be engineered via optical 
pumping. 

A. Transverse-fleld Ising model (TFIM) with 
spontaneous emission in the eigenbasis of the field 

In this subsection, we consider a driven-dissipative 
model described by the Hamiltonian (written in terms 
of Pauli matrices (t“) 

= + (13) 

(b) » 

where the first term with J > 0 [104] is the nearest- 
neighbor interaction on a d-dimensional cubic lattice with 
d = 2 or 3, and with the dissipation via the Lindblad op¬ 
erator at each site = -x/To-“ = -x/T {af — which 

is simply the spontaneous emission from spin up jf) to 
spin down ]),). In the absence of dissipation, the ground 


state of Eq. (13) is known to give rise to a quantum phase 
transition [ 1 ] from an ordered phase (completely ordered 
at J/A ^ 1), where the Z 2 symmetry ax —>• —ax is 
broken and {ax) ^ 0 , to a disordered phase (fully disor¬ 
dered at J/A <C 1), where the symmetry is restored and 
(ctx) = 0. Eor the dissipative system considered here, Z 2 
is still a symmetry of the master equation. To see this, 
consider the transformation 

cTx -t -CTx, a-y -o-y, a^ ax, (14) 

which respects the noncommutative algebra of Pauli 
operators, and under which Li —>■ —Li. Since the 
master equation is bilinear in Li and L|, the overall sign 
(or phase) of the Lindblad operators is insignificant. 
One may then expect a phase transition between an 
ordered and a disordered phase even in the presence of 
dissipation. We will first briefly discuss the mean field 
(ME) prediction, and then compare and contrast it with 
the more careful Keldysh field-theoretic treatment. 

Mean field —Eor (3 J — A)A > 0 (with 3 = 2(i the 
coordination number), as P is increased, the ME gives 
a continuous transition from the ordered phase (erf) = 
const ^ 0 to a disordered phase (erf) = 0 (see App. 
A). Most qualitative features of the ME are confirmed 
below by the field-theoretic treatment with quantitative 
corrections. One qualitative difference, however, is that 
the second-order transition seems to be replaced by a 
first-order transition for sufficiently strong dissipation. 
Eurthermore, the Keldysh formalism reveals that the 
phase transition belongs to the universality class of 
the classical Ising model with an effective temperature 
determined by the microscopic parameters in the model. 
Interestingly, the effective temperature remains finite 
even in the limit P —^ 0 . 

Field theory —For a proper field-theoretic treatment, 
we first write the spin operators in terms of hard-core 
bosons 

(^i=ai, cr+ = al, a- = 2ajai - 1. (15) 

The hard-core constraint can be implemented via a large 
on-site potential U ajajaiai with U —>■ -l-oo. We are 
particularly interested in taking the continuum limit via 
tti a(x), where the operator a(x) varies continuously 
in space. Since the ME predicts uniform phases, our as¬ 
sumption should be justified as a starting point for the 
field theory. However, as is typical with the transition 
to the continuum, hard-core features become soft in the 
continuum after coarse graining where short wavelength 
modes are eliminated. For example, the (near-)critical 
behavior of the classical Ising model with discrete values 
Si = ±1 is mapped to the field theory [ 1 , 82]. At 
long wavelengths, the Ising spins can be coarse-grained 
to a continuous field cj), and interact via a soft term 
which respects the original symmetry of the problem 
[erf —>■ —af reflected by (/(x) —>• —(/(x)]; in principle, 
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one must also include higher-order terms ( 0 ®, etc.) that 
respect the symmetry, but such terms are less relevant 
under RG and can be dropped. The quantum TFIM also 
enjoys a similar mapping to the (f)^ theory in the con¬ 
tinuum albeit in one higher dimension [1]. With this in 
mind, we shall frequently use the mapping to the con¬ 
tinuum, and add the quartic term introduced above with 
a finite strength of the interaction. The long-distance 
behavior of our model should be insensitive to this as¬ 
sumption. This is especially the case near the critical 
point where (a) and (a^) are small, and a phenomenolog¬ 
ical expansion and truncation of the interaction term— 
consistent with the underlying symmetries—at the quar¬ 
tic order is further justified. 

In the continuum, the first term in the Hamiltonian 
(13) written in terms of the bosonic operators should be 
expanded up to the second derivative in spatial coordi¬ 
nates; higher derivatives can be ignored in the long wave¬ 
length limit. Other terms in the Hamiltonian map to the 
continuum in a straightforward way. The full Hamilto¬ 
nian then reads 

H = J -J [a(x) -I- a'l'(x)] -I- [a(x) -|- a'''(x)] 

-I- 2Aa^(x)a(x) -|- [/a^(x)a^(x)a(x)a(x), (16) 


with the Lindblad operator Lx = ■\/ra(x); the lattice 
spacing is taken to be unity for simplicity. As we shall 
see, this continuum model reproduces the MF equations 
for small dissipation (cf. App. A) but indicates the failure 
of the mean field for large dissipation. 

Next we map the master equation to the Keldysh 
path integral in terms of classical and quantum fields 
aci,q(t,x). It is more convenient to work with the real 
and imaginary parts of Ud/q = Cd/q + idd/q- Equations 
(5, 6 ) then yield the Keldysh Lagrangian (cf. Ref. [44]) 


approximation as [105] 

4r p ~ ~ F^ 

^eff “ ~^CqOtCcl 4t/ CgV CqI (</ ^ ^CqCd 

+ iT{l +— —-^)c^dC-q (18) 

+ 0{cqdtCd, i{dtCqf, Cdcl) , 

where 0{-) contains irrelevant terms in the sense of RG 
due to their higher-order time derivatives, or higher pow¬ 
ers of the quantum field (see Sec. H). We briefly remark 
that, in the absence of dissipation, the linear time deriva¬ 
tive vanishes, and one should keep the second-order time 
derivative, which, combined with the second-order spa¬ 
tial derivatives, gives rise to the field theory in one 
higher dimension , and the critical point corresponds to 
setting the mass term (the coefficient of CdCq with F — 0 ) 
to zero. In the presence of dissipation, the linear time 
derivative is more relevant, and the quantum criticality 
is inaccessible. In this case, the critical point obtained by 
setting the mass term to zero matches exactly with the 
MF equation derived in App. A. However, we will see 
shortly that there are important caveats indicating the 
failure of the mean field, especially for sufficiently strong 
dissipation. 

Neglecting the irrelevant terms, the quantum field Cq 
appears at most quadratically in the Keldysh action, and 
thus an exact classical Langevin equation emerges for Cd 
{cd -t c): 


4r 

- -^9tc = [-4J + r + uc^{t,x.)] c{t,x) + ^(t,x). 


with the parameters 


r = — J + A 


A 


u = U\l- 


aW ’ 


(19) 

( 20 ) 


where the noise ^ is correlated as [57, 61] 


/Ik — 2 ( Ccldtdq ~\~ dcldtCq) 4(/ Cq\/ CqI J CqCd 
^(CgCc/ dddq') T^i^Cddq Cqdd^ ^q) 

~U{cd + Cq -b d?d + dg){cdCq + dddq), (17) 


where J = A = 4A, and 3 = 2(i is the coordina¬ 
tion number. Notice that there is no gradient term in 
the fields dd/q, which makes them purely local in spa¬ 
tial coordinates. Furthermore, they are “gapped” due to 
the terms — Adcdq-fiFdq in the action (while Cd/q can be 
tuned near criticality), but we cannot just drop the dd/q- 
dependence as there would be no time derivative acting 
on Cd/q- We thus integrate out the former, and find an 
effective action for the field Cd/q- Since dd/q are gapped, 
they can be integrated out at the level of the quadratic 
terms in the first and second lines of Eq. (17); the ef¬ 
fective Lagrangian simply follows from the saddle-point 


(e(t,x)C(t',x')) = 2 r(^l + ^]) 5(x-x')<5(t-t'). (21) 

Equation (19) describes dynamics that is mathematically 
equivalent to the dynamics near the thermodynamic equi¬ 
librium of the field c. In this sense, the steady-state so¬ 
lution to Eq. (19) is simply given by the thermal Gibbs 
ensemble (normalizing u) 


P[c(x)] ^ exp 




2 


( 22 ) 


with 


Teff = 


A 

T 



(23) 


where Teff is obtained by imposing the fluctuation- 
dissipation relation discussed following Eq. (12). We 
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point out that the effective temperature goes to a fi¬ 
nite value even for F —>■ 0. Therefore, even with in¬ 
finitesimal dissipation, the system, contrary to what one 
might naively expect, does not get arbitrarily close to the 
ground state of the effective Hamiltonian in Eq. (22). At 
the same time, for infinitesimal dissipation, the approach 
to the steady state will typically take a very long time. 

The partition function defined as the sum over all con¬ 
figurations of c weighted by the probability distribution 
(22), f D[c] P[c], is nothing but the 4>‘^ theory, which de¬ 
scribes the universality class of the classical Ising model. 
The critical point is given by r = 0, at which point we 
have J = 43 J = A(l-|- F^/A^). Comparing against 
Eq. (23), it becomes evident that, at the critical point, 
J/TeS = 1/3- Interestingly, the same relation also de¬ 
scribes the mean-field solution for the critical point of 
the classical Ising model H = —JJ2{ij) without the 
transverse field and dissipation. This implies that the 
transverse field together with dissipation play the role of 
thermal fluctuations, where even the value of the critical 
temperature is matched (at, and most likely also beyond, 
the mean field level) with the classical Ising model. This 
surprising feature is most likely related to the fact that 
both the transverse held and dissipation are dehned in 
the cr^ basis. Comparing to Eq. (13), one can see that 
the effective result of dissipation is simply to reduce the 
original Hamiltonian to one with only the Ising term, 
which, however, should be regarded at a hnite effective 
temperature. On the other hand, in all other models 
that we study in this paper, the hnal effective Hamilto¬ 
nian and the corresponding thermodynamic model bear 
no such obvious relationship to the original dissipative- 
driven model. 

For our held-theoretic treatment to be valid, the par¬ 
tition function should be convergent, and specihcally 
M > 0, i.e. F < 4A. For F > 4A, the sign of the quartic 
term is negative, and the model described by Eq. (22) 
exhibits an instability towards a phase where c —^ ± 00 . 
Of course, the sum (trace) over spins is always conver¬ 
gent (spin excitations will be saturated), and thus there 
should be higher order terms such as c® with a positive 
coefficient in Eq. (22) making the functional-integral fully 
convergent. In this case, one finds that the second-order 
phase transition is replaced by a first-order phase tran¬ 
sition in a model described by the effective Hamiltonian 
density H ~ 2J(Vc)^ -I- {rl2)(? + uc'^ + nc® where u < 0 
and n > 0, see p. 173 of [83]. 

Dynamics: The Langevin equation (19) clearly indi¬ 
cates that the dynamics is diffusive, i.e. the dynamic 
exponent is ^ = 2 at the mean-field level with correc¬ 
tions due to fluctuations captured by loop diagrams. In 
short, we have reduced the starting dissipative spin sys¬ 
tem to the dynamical Landau-Ginzburg model where the 
dynamical field is not conserved. The latter model falls 
under the so-called model A dynamics of the Halperin- 
Hohenberg classification [84], where 2 ; is known via 
epsilon-expansion or numerical evaluation. 


B. Isotropic XY model with coherent drive and 
spontaneons emission 

In this subsection, we consider the lattice Hamiltonian 

H = - + h.c. -f H ^ erf -f A ^ uf , (24) 

(y> * * 


assuming nearest-neighbor interaction with J > 0 [106] 
on a d-dimensional cubic lattice with d = 2 or 3, together 
with the dissipative process via the Lindblad operator 
at each site Li = vT a~. Without the coherent drive, 
H = 0 , the dissipation drives the system to a “dark” state 
where all spins are in state j),); this happens because 
this dark state is an eigenstate of the Hamiltonian for 
H = 0. To find a nontrivial steady state, we consider a 
finite coherent drive VL ^ Q. 

Mean field —The phase diagram of this model, at the 
level of ME, includes a region with a uniform ((erf) = 
const ^ 0 ) stable phase, and a bistable regime, in which 
two uniform stable phases exist (see App. B). The ME 
predicts a critical point at which a continuous transition 
occurs between stable and bistable regions. While the 
full ME phase diagram has more structure to it, we 
are particularly interested in this criticality and closely 
investigate its vicinity. We show that the continuous 
phase transition is similar to that of the Ising-type 
liquid-gas phase transition, but find that bistability is 
an artifact of the ME. 

Field theory —We begin by writing the spin operators 
in terms of hard-core bosons, Eq. (15), and implement the 
hard-core constraint via an on-site quartic potential as in 
the previous subsection. We remark that the parameters 
in the model can be chosen such that (a) and (al) are 
small near the critical point, which further justifies the 
truncation of the interaction at the quartic order. In the 
anticipation of a uniform phase, we cast the Hamiltonian 
in the continuum as 

H= f — Ja^(x)V^a(x)-I-A a^(x)a(x) 

J X 

-I- H (a(x) -I- a^(x)) -|- U a^(x)a^(x)a(x)a(x), (25) 

where A = 2A — 3 J with 3 the coordination number; we 
shall drop the tilde for notational convenience, but always 
mean A in the rest of this subsection. Similarly, the 
Lindblad operators go to Li —>■ •\/ra(x). The Keldysh 
action is then given by 


Sk= 


0 


idt+JV^-A-f Wuci 


^zdt-bJV^-A-bf 
U 


iV 


- V2il{aq + a*)- — (jaczl^ + [a,]^) (acia* + c.c.) . (26) 


Before going further, we compare the continuum de¬ 
scription in Eqs. (25,26) with the original lattice model. 


The two descriptions yield almost identical MF equations 
with a similar critical behavior and exhibit stable and 
bistable regions (MF in the continuum will be discussed 
shortly); however, for the corresponding regions to match 
exactly, U is to be substituted in terms of the original pa¬ 
rameters in the lattice model. For the sake of clarity and 
to avoid confusion, we shall regard the continuum de¬ 
scription in Eqs. (25,26) as our fundamental model and 
a starting point for further investigation. In fact, closely 
related models arise naturally in the context of interact¬ 
ing Rydberg polaritons in free space [85, 86 ]. However, 
we maintain that, at least for the regions in the param¬ 
eter space with small excitation density, the continuum 
description should be a good approximation to the lattice 
model. 

We first embark on a detailed study of the vicinity of 
the MF critical point. To this end, the mean field solution 
is first derived for the continuum model via 8Sk/ 8a* = 0 , 
yielding 


(^-A + z0ao-V^H-ic/|aoPao = 0, (27) 

where uq = (ad) is the MF value of the classical field. 
Depending on the parameters, this equation has either 
one stable solution, or three solutions, only two of which 

I 


are stable. Near the critical point, these solutions are 
continuously connected. To characterize this point, we 
first define ( = —A/U, 7 = F/C/, o = H/{7, and n = 
|aop /2 the excitation density (the factor of 1/2 appears 
due to the definition of ad in terms of the original fields). 
Equation (27) can then be cast as 


(c-.7 + (|)‘ 


(28) 


At the critical point, the above parameters are given by 




(29) 


which is easily verified by noting that the three roots of 
Eq. (27) become degenerate and give the critical excita¬ 
tion density of 



(30) 


As the next step, we expand the Keldysh action, Eq. (26), 
around the MF solution in Eq. (27), i.e. ad —>"00 + o,ch 
and find 


r I 0 D#^ 2 (^,k)\ 

/ ^'''(w,k) A(a;,k) 

Wx2(^,k) ) 

I (2ao|flci|^ + Uq ®c;) A c.c. -I- (jaczj^ -I- |ng|^) {ada^ c.c.) 

J t.X 


Sk = 

u 


(31) 


where A(w,k) = ^aci(w,k) a*;(— w,— k) aq(oj,k) a*(—oj,—k)J , and 

^^x 2 (w,k) = \ and 


^ 2 ^'^ are 2 x 2 matrices given by 0 ^x 2 = 


iT 0 ’ 
. 0 fF, 


^2x2 (^) k) — 


'(jj + *F/2 — Jk^ — A — U\aoY 


-{U/2)a*o 


*2 


-{U/2)al 

-to - iT/2 - Jk^ - A - U\ao\ 


(32) 


The second line of Eq. (31) includes cubic and quartic 
terms in the action. To find the dissipative spectrum of 
fluctuations, or, more precisely, the relaxation rate, we 
solve det k) = 0 and find 

Wk = -^^± V(A + C/|ao|2 + Jk2)2 - (t//2)2|ao|4 . (33) 

It is easy to see that one of the two eigenvalues vanishes at 
the critical point as expected. To approach the critical 


point, we can fix n = (or equivalently fix ag) and 
o = 0 ( according to Eqs. (29,30), and take 7 —>■ For 
7 > 7 ^, there is a unique solution to Eqs. (27,28), while, 
for 7 < 7 f, two stable solutions continuously emerge. 
Casting uj in units of U, and working in the long wave¬ 
length limit, we find the relaxation rates as 
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where the coefficient of in the first eigenvalue is pro¬ 
portional to J but is normalized to unity (by rescaling 
space), while, in the second eigenvalue, the momentum 
dependence is entirely neglected due to the dissipative 
gap. We are further interested in finding a natural basis 
for ^^x2( w, k) in order to break the fluctuations of the 
field into massless and massive components. To this end, 
we first drop the k-dependence in Eq. (32) as we are in¬ 
terested in the long-wavelength limit; we will deal with 
the k^ term later. (On similar grounds, we can also drop 
w, but it does not further simplify our task.) Also, to 
simplify computations, we absorb the phase of ag in the 
definition of thus ao,ag |ao| in the off-diagonal 

elements of the matrix in Eq. (32) as well as the coeffi¬ 
cient of the cubic term in Eq. (31). The resulting matrix 
^ 2 x 2 takes the form (in units of U) 




^w-kf7/2-C/3 -2C/3 

^ -2C/3 -w-17/2-C/3 


(35) 


cast the interaction terms in Eq. (31) in the new basis. 
We start with the cubic interaction, 

A [chdq + dhcq] . 

The precise value of the coefficient A will not be impor¬ 
tant. Next, the quartic term expanded in terms of Cd/q 
and dci/q generates many terms (with it > 0): 

(Cci + + Ccldcl — Cqdq) 

^ dcldq -f ‘2Cqdcl ‘^Cddq^ . 

And finally the noise term (in units of U) takes a simple 
form in the new basis iyjV'gP ~ *7 {cq + d^ — Cqdq). In 
writing the full Keldsyh action, we note that the scal¬ 
ing dimension of the fields renders terms with two or 
more powers of the quantum field Cq irrelevant (except 
the noise term), as explained in Sec. II. We keep the rel¬ 
evant terms in the field Cd/q, but also include some (not 
all) cross terms with dd/q for reasons that will become 
clear shortly. The Keldysh action then reads 


To simplify the form of the action, we change the basis 
to 


Cd/q{^, k) = =F e^'-^^^GdjqiuJ, k) -H 
dd/qM = e^"/^,7,(w,k) + 


ci/g(-‘^>-k)J , 

(-w,-k), 

(36) 


which allows us to rewrite the quadratic part of the 
Keldysh action as (for the moment, not including the 
k dependence in the kernel and the quadratic quantum- 
noise term) 


C 


( 2 ) 

K 


^ Cq{-U},-k)Cd{uJ,k) 

+ dq{-uj, -k)dd{_uj,k) 


7 , C ' 

2 + 7i, 

2 V3 


(37) 


where we have neglected a multiplicative constant. Equa¬ 
tion (37) clearly indicates that Cd/q can be tuned near 
criticality by taking 7 —^ 7c, while dd/q are always 
gapped. Also notice that it follows from Eq. (36) that 
the fields Cd/q{t,^) and dd/q(t,x.) are real-valued. In¬ 
verting the basis in Eq. (36), and casting it in real space 
and time, we find 


ad/qit,^) = (Te'’'“/^Cci/g(t,x) -h x)) . 

(38) 


With this representation, we can now cast various terms 
in the action in terms of the new fields. The gradient 
term a*V^ac; -I- c.c., i.e. the k^ term in the action, now 
takes the form (with a normalized coefficient) 

CqV'^Cd H-, 

where the ellipses denote dci/g"dependent gradient terms, 
which are ignored since dd/q are gapped [107]. Next we 


Ck ~ Cq{-dt - r + V'^)Cd - U CdCq + *7 Cg 
+ dq{-dt - R)dd +'ild^q 
“f (cciC^g T c^cfCg) — 2u Cqd^i 3u c^idddq -!-■■■ 

- i 7 Cqdq , (39) 


where 


7 C 
2 73’ 


^-2 +VS' 


(40) 


We have organized Eq. (39) into (the first line) terms de¬ 
pending on Cd/q only, (the second line) quadratic terms 
in dd/q, and (the third and fourth lines) various dd/q- 
dependent nonlinear and cross terms; the ellipses denote 
nonlinear terms not written explicitly. Also we have not 
fully kept track of various coefficients (except those of r 
and R) since they will be inconsequential for our conclu¬ 
sions. 

If we simply drop the dcz/g-dependent terms, we get 
for the Keldysh action the first line of Eq. (39). By tak¬ 
ing the second step of writing the corresponding Langevin 
equation and subsequently mapping to the partition func¬ 
tion as outlined in Sec. II (see also the previous sub¬ 
section), we would find the thermodynamic Landau- 
Ginzburg theory with the Z2 symmetry, and the as¬ 
sociated second-order phase transition as r —>■ 0. Taking 
into account the fields dd/q and their fluctuations, we 
next show that the Z2 symmetry is spoiled, but a contin¬ 
uous phase transition emerges akin to that of the liquid- 
gas transition. To start with, let us consider the effect of 
fluctuations to first order in A; at this order, we find 


^{dd)cq, 


(41) 


generated by integrating out dd in the cubic term. This 
term acts like a “magnetic” field, and breaks the Z2 sym¬ 
metry of Cd,q —t —Cd,q- However, it can be absorbed into 
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the parameter Q in the original model. In fact, fluctu¬ 
ations can modify the position of the critical point, and 
thus the above equation can be regarded as the correc¬ 
tion to the MF position of the phase transition. We thus 
consider the effect of fluctuations at higher orders in A 
and u. The nonlinear terms in (the third line of) Eq. (39) 
expanded to second order generate various terms; of par¬ 
ticular interest to us is the term proportional to 


Xu 


Xu 


^ci^q 


^ci 


{dliit)dq{t + t)) 

-3 {dci{t)dq{t)dli{t + t)) I 
/'drG«(r) [G^(0)-G^(r)] , (42) 


with all the fields evaluated at the same spatial co¬ 
ordinates since there is no gradient term in d^i/q- 
Also G^{t) = —i {dci(t)dq{t + t)) and G^(t) = 

—i {dci{t)dci{t + t)) are the retarded and Keldysh 
Green’s functions, respectively, for the fields dcZ/iji which 
in Fourier space become G^(uj) ^ {iuj — R)~^ and 
iG^{ijj) ~ G'^(u;) 7 G^(w)*. These functions have a 
nonzero support only for r < \/R ^ 7~^, which is why 
the vertex in c^i/q in Eq. (42) is approximated to be local 
in time. Due to the nonvanishing integral in this equa¬ 
tion, fluctuations will generate a term of the form ^ chcq 
in the action with some constant v. With the first- and 
cubic-order terms generated as the result of fluctuations, 
there is no apparent Z2 symmetry; the full partition func¬ 
tion then takes the form {cd —>■ c) 


1 

7 

(43) 

where the exact coefficients of various terms are disre¬ 
garded [108]. Similar steps to those outlined in Sec. II 
(see also Sec. Ill A) are taken to obtain the partition func¬ 
tion from the Keldysh action and the resulting Langevin 
equation. Terms with odd powers of c should be traced 
back to the fact that the full action (31) as a function(al) 
of {cci,5, dcZ,g} is not symmetric under the simultaneous 
transformations Cd^q —>• —Cd,q and dd,q —>■ —dd,q- DO" 
spite this fact, we can absorb the linear term in the pa¬ 
rameter D, and shift the field c by a constant (c —>■ cq -be) 
to eliminate the third-order term; the constant term 
merely modifies the MF critical point around which we 
have expanded the action. A similar scenario arises in 
the liquid-gas phase transition, where there is no obvi¬ 
ous symmetry, however, one can choose parameters such 
as density to eliminate odd terms. This phase transition, 
despite the absence of symmetry, is of the Ising type [83] . 
We thus conclude that the driven-dissipative model con¬ 
sidered in this subsection undergoes a continuous Ising- 
type phase transition. We further remark that there is 
no true thermodynamic bistability; the true steady state 
of the system is given by the minimum of the exponent 
in the partition function (43), which is unique for generic 
values of h and v. We stress that our argument does not 


/ (Vc)^ + hc + rc^ + vc^ + 1 

J X 



rely on a similar minimum criterion in thermodynamics, 
but is simply derived in our nonequilibrium model; in the 
thermodynamic limit, a minimum of the exponent [in Eq. 
(43)] , being extensive in the system size, is infinitely more 
likely to occur than any other state including other local 
minima of the exponent. 

Similar models describing a dissipative gas of Rydberg 
atoms have been studied in Refs. [28, 31]. A notable dif¬ 
ference is that the interaction in Refs. [28, 31] is a long- 
range interaction of the Ising type (in contrast to 

the nearest-neighbor flip-flop interaction in the present 
subsection). Nevertheless, the corresponding MF analy¬ 
sis performed in Ref. [31] is almost identical to the MF 
equation of this subsection. Therefore, at least in uni¬ 
form phases, the model in Ref. [31] might be amenable 
to a similar treatment. It would be interesting to study 
the effect of fluctuations beyond the mean field; see also 
the comparison with an approach based on a variational 
principle for steady states [32]. Finally, we note that our 
considerations here should be directly applicable to uni¬ 
form phases of driven-dissipative Bose-Hubbard models 
and their critical behavior studied in Refs. [30, 87]. 


C. Anisotropic XY model with spontaneous 
emission 


In this subsection, we consider the Hamiltonian 




2d 




<b> 


(b> 


, (44) 


assuming nearest-neighbor interactions with J > 0 on 
a d-dimensional cubic lattice with d = 2 or 3, together 
with dissipation via the Lindblad operator at each 
site Li = -x/r a~. While spontaneous emission tends 
to create a state where all spins point down, such a 
state is not an eigenstate of the Hamiltonian. The 
interplay of the dissipation and the effective drive in the 
Hamiltonian can give rise to nontrivial steady states. An 
experimental realization of this model using ultracold 
atoms in the ground electronic state weakly dressed with 
highly excited Rydberg states is proposed in Ref. [29]. 


Mean field —The MF is studied in Ref. [29] for the 
more general XYZ model with ^ Jy ^ J^. For 
sufficiently weak spontaneous emission in the model 
considered here, the MF (see App. C) predicts a sponta¬ 
neous symmetry breaking that gives rise to a staggered 
XY steady state with spins on neighboring sites pointing 
in different directions. For larger values of spontaneous 
emission, one finds a disordered paramagnetic state. 
Our field-theoretic treatment confirms this picture in 
three dimensions, however, in two dimensions, we show 
that the XY phase cannot be realized. This happens 
because the effective temperature, emerging due to the 
dissipation, is larger that the Kosterlitz-Thouless tem¬ 
perature associated with the transition from short-range 
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to algebraic long-range order in two dimensions. 

Field theory—As in the previous subsections, we first 
map spins to hard-core bosons, and represent their hard¬ 
core nature via a quartic term. This mapping is par¬ 
ticularly good close to the phase transition between the 
paramagnetic and the XY phase because cr^ « — 1 near 
the phase boundary. In anticipation of the staggered XY 
phase, we recast the Hamiltonian on two checkerboard 
sublattices A and B as 

H = ^^(&la]-k6»aj)-kC/ + U ^ blbjbA , 

(ij) iSA jeB 

(45) 


where J = 2J/d. Note that we have chosen the same 
coefficient for the interaction terms on both sublattices 
to make manifest the symmetry of the underlying spin 
model. The Keldysh Lagrangian at the quadratic level, 
and in Fourier space, can be written as (the integral over 
frequency is limited to w > 0) 


J 






+ J(k) [a*{-u;, -k)6*;(w, k) -h aci{-uj, -k)6,(u;, k) -h c.c.] , 




(46) 


where we have defined >^(k) = 

J [cos(fca,) -I- cos(fcy) -I -], with dots standing for 

cos{kz) in the case of d = 3 dimensions. A quadratic 
Keldysh action with the assumption that spins map to 
soft-core bosons, and an ansatz of a uniform phase (i.e., 
by identifying a and b fields), has been first constructed 
in Ref. [29]. The above quadratic action becomes 
diagonal in the basis defined by 


basis, we invert the eigenbasis in Eq. (47) and cast it in 
space-time coordinates, 


, , 1 

Te^*"/^c:,/,(t,x)+e±-/V:,/,(t,x)' 

bd/q{t.^) = ^ 

Te'^""''''^Cd/g(t, x) -k e^"^^^dd/q{t, x) 


(49) 


Cci/qi^, k) = =F k) -f i, -k) 


(47) 


which casts the quadratic part of the Keldsyh Lagrangian 
into 

=\{ ^ + d(k)]c*(a;,k)cci(a;,k)-hc.c. 

+ [iu)-^- J(k)]d*(w,k)dci(w,k) 4- c.c. 

-k iT (|cg(w,k)p -4 |d,(a;,k)p) j. (48) 

At long wavelengths, k —> 0, and J(k) —>• J(0) = 2J > 0. 
Therefore, the fields c^i/q can be tuned near criticality 
with a vanishing gap characterized byr = r/2— 2J, 
while dci/q are massive with (the imaginary part of) the 
gap i? = r/2 -I- 2J, and can be integrated out; however, 
we must not drop them before considering the (nonlin¬ 
ear) interaction terms. Casting the quartic interaction 
terms in the continuum and then in the Keldysh basis, we 
get (-t//2) /j ^ (Iflcip -k |a,p) [acia* + c.c.) plus a simi¬ 
lar term for 6^;/^. To write the interaction in the new 


The interaction then takes the form (up to a multiplica¬ 
tive constant) 

U f (|cc;P + jdcZp + |cqrp + |dgp) (^—Ccid* + ddC* + c.c.) 

J £,x 

- {Ccidh + Cqd* - c.c.) {ccic* + dcid* - c.c.) . 

(50) 

If we simply drop dd/q, there will not be any nonlinear 
terms. Instead we should find the new interaction ver¬ 
tices generated via integrating out dd/q- We thus expand 
the Keldysh functional integral to the first few orders in 
the interaction. To first order, the effective interaction 
vanishes once averaged by the Gaussian functional inte¬ 
gral due to Eq. (48) since all the terms in Eq. (50) are 
odd in dd/q- To second order, we generate the vertex (up 
to a positive prefactor [109]) 




\^cl\ ^cl^q 


iG^ {0) I dr G^{t) . 


- t,X 


(51) 


Many terms contribute to this vertex, but its precise coef¬ 
ficient is immaterial for our considerations. As in the pre¬ 
vious subsection, the fields are evaluated at the same spa¬ 
tial coordinates (the gradient term in dd jq is ignored due 
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to the dissipative gap). Also the retarded and Keldysh 
Green’s functions for the fields d ci/q are, similar 

to the expressions given before, G^{u)) ~ [ioj — and 
iG^(uj) ~ G'^(w) r G^(a;) with a nonzero support for 
T < r“^, hence the local form of Eq. (51) in time t. 
In short, fluctuations will generate a term of the form 
—u\cci\^Ccic* in the action; the facts that J drG^ij) = 

G^{uj = 0)- l/R < 0 and iG^ {t = 0) - T/R > 0 

ensure that u > 0. Various other terms generated by in¬ 
tegrating out d^i/q either give corrections to the existing 
terms in the action, or produce irrelevant terms in the 
sense of RG. Importantly, the effective action respects 
the 17(1) symmetry (c —>■ ce*®) beyond the quadratic or¬ 
der. The final form of the Keldysh Lagrangian, with the 
relevant terms only, becomes 

C-K =^|c*(-i9t + ^JV^-r)ccz-f c.c. 

- u\Ccl\^ [cclCq + c.c.) -b irlcql^j. (52) 

The quantum vertex appears at most quadratically, lead¬ 
ing to a classical Langevin equation with a noise term 
which can be interpreted as an effective temperature. 
The corresponding steady state is then described by the 
thermodynamic partition function (ccz —>■ '0) 


ll[0(x)] exp 


1 

Tea 


L 


J |V0p + r 101^ + 


(53) 

where Teff = T -b • • • with the ellipses being the cor¬ 
rections due to renormalization; the effective tempera¬ 
ture and its precise coefficient (unity) are obtained using 
the fluctuation-dissipation condition. The partition func¬ 
tion (53) belongs to the universality class of the classical 
XY model. One should then expect off-diagonal order in 
d = 3 dimensions, and a Kosterlitz-Thouless (KT) tran¬ 
sition in d = 2 dimensions. Nevertheless, we remark that 
a possible emergence of the Kardar-Parisi-Zhang equa¬ 
tion [88] similar to Ref. [45] may invalidate the above 
analysis in two dimensions; however, we show that, even 
in the absence of such mechanism, the constraints on the 
Kosterlitz-Thouless temperature would prevent the sys¬ 
tem from realizing the XY phase in two dimensions. To 
see this, we should examine the condition for the KT 
transition. Denoting 0 = j0je®® in the ordered phase, 
the algebraic long-range order is realized when (J = J in 
d = 2 dimensions) [61] 


^ 2^ 

Teff TT 


(54) 


Normally, at sufficiently low temperatures, this condition 
is well satisfied; however, in this case, as Tea 0, or 
equivalently T —>■ 0, we also find—from the mean-field 
analysis—that 0 —>■ 0 in the same limit. To see this, note 
that 

0 ^ c = e^/^a -b ^ ^{a^ + a^). (55) 

v2 


Therefore, j0j^ = l(c)j^ = where we have taken 

(a^) = (u^); the value of (a^) is inserted from the mean- 
field analysis in two dimensions. 


{an 


04 JT - r2 
4J 


(56) 


With these expressions, and Tea ~ E, the condition (54) 
takes the form 


1 1 1 

4 16j ^ TT ’ 


(57) 


with j = J/r being larger than 1/4 in the XY phase, a 
condition that is never satisfied, and becomes even worse 
with increasing J (or decreasing P). Of course, in evalu¬ 
ating the above expressions, we have used mean-field ex¬ 
pressions which can be modified, and relied on our field 
theory description away from phase boundaries. How¬ 
ever, one should expect that the algebraic long-range or¬ 
der in two dimensions will be significantly diminished, if 
not completely disappear. 

Dynamics'. In d = 3 dimensions, the Langevin equa¬ 
tion corresponding to Eq. (53) indicates that the dynam¬ 
ics is diffusive. The dynamical field is not conserved, and 
thus the dynamics falls under Model A of the Hohenberg- 
Halperin classification for the Landau-Ginzburg model 
with N = 2 components [84]. 


D. Isotropic XY model with incoherent pumping 
and interaction-induced loss 

In this subsection, we consider the Hamiltonian 

id = - J ^ cr+ a~ + h.c. + A ^ cr/ , (58) 

(ij) » 


with J > 0 on a three-dimensional cubic lattice [110]. 
In all the models in the previous subsections, we have 
chosen a simple dissipative process where spins at 
different sites spontaneously and independently decay 
from jt) to jj,). In this subsection, in addition to 
the spontaneous emission via L\ = vT a ~, we also 
consider pumping defined via = -^/Ppcr/", and an 
interaction-induced loss described by the Lindblad 
operator Lfj = \/^a^a~a~ for nearest neighbors i and 
j. The latter is an example of a more complicated type 
of dissipation that depends on the correlation between 
nearby sites. In this case, the jump operator checks 
if there are two excitations on neighboring sites, and, if 
there are, kills one. This type of operator is natural in 
systems where a particular laser coupling scheme creates 
dark states, or pseudo-spin states, as linear combinations 
of the microscopic energy levels; however, interaction 
between neighboring pseudo-spin states shifts them out 
of resonance, and can lead to the decay of one of them 
[18, 89, 90]. 







13 


Mean field —Without the interaction-based loss, 
this model is rather trivial. For a fixed decay rate F, 
the excitation density increases with increasing Fp, but 
(cr“) = (cr^) = 0 since the steady state can be easily seen 
to be a product of single-site density matrices diagonal in 
the az basis. Specifically, for Fp = F, the system is in an 
infinite-temperature state. At the level of the MF, the 
same qualitative behavior persists even in the presence 
of the interaction-induced loss, although it changes 
quantitatively. More importantly, the phase predicted 
by the MF does not break the U{1) symmetry of the 
problem, i.e. {a^) = (a^) = 0; see App. D. However, we 
shall see that the Keldysh path-integral approach gives 


a continuous transition from the disordered phase to a 
phase with spontaneous continuous symmetry breaking. 

Field theory —We start by assuming that k ^ F, Fp. 
With this assumption, the excitation density is rather 
small, and one can safely represent the spins in terms of 
soft-core bosons but with quartic interactions as in the 
previous subsections. We then find the Lindblad oper¬ 
ator Ljj = OjUjai- In the continuum description, this 
operator can be cast as —>• a’l’(x)a^(x) plus gradient 

terms which are less relevant [111]. Some algebra yields 
the Keldysh Lagrangian (with a normalized J) 


C-K = 




\idt + i — 2 ^ 


idt + + pL — i 


*(r + rp) 


• r-r 



-f C.C.^ C.C.) , 


(59) 


where we have introduced p. using the freedom (due to 
the symmetry) in choosing a rotating frame with (a) ~ 
exp{—iu}Qt). In writing the Lagrangian, we have ignored 
terms at the quadratic or higher orders in the quantum 
field—except the noise term ijogp—as they are irrelevant 
in the sense of RG. Casting the Keldysh path integral as 
a Langevin equation, we find (ad ip) 


idt + + p + 

= f(^,x), 


{7|V'P + '0(i,x) 


(60) 


with 


m x)r (t', X')) = (F + Tp)dit - t')Six - x'). (61) 

Before considering fluctuations, we study the mean field 
solution at the level of the Keldysh action or the corre¬ 
sponding Langevin equation (60). We shall see that even 
the mean field, at the level of the path integral, improves 
upon the MF on the lattice model. In the absence of fluc¬ 
tuations, a uniform phase exists provided that Fp > F, 


limit. Our model is slightly different because the dis¬ 
sipative nonlinearity arises at the fifth, rather than the 
third, order in Eq. (60); however, the RG procedure (and, 
most intuitively, momentum-shell RG) creates all pos¬ 
sible terms consistent with symmetry. Therefore, our 
model flows to the same universality class as the model 
in Refs. [43, 44], or, equivalently, that of Eq. (53) with 
r = (F — rp)/2 and Tgg = F -|- Fp. Specifically, the gradi¬ 
ent term, being purely coherent at the level of the original 
Hamiltonian, becomes dissipative in the course of RG. 

Similar considerations apply to a different model de¬ 
scribed by the Lindblad operator Ljj ^ '^7 which 

kills both excitations on neighboring sites. The lattice 
MF analysis again fails to capture the full phase dia¬ 
gram. The field-theoretic treatment in this case becomes 
almost identical to the model in Refs. [43, 44] which then 
predicts a phase with spontaneous continuous symmetry 
breaking. 


IV. CONCLUSIONS AND OUTLOOK 


iP = 



(62) 


for a constant phase 0; the real part of the bracket in 
Eq. (60) vanishes by appropriately choosing the con¬ 
stant p. The solution in Eq. (62) explicitly breaks the 
17(1) symmetry. In fact, a similar model was studied in 
Refs. [43, 44] where the authors concluded that the sys¬ 
tem becomes purely dissipative under RG, and specif¬ 
ically the real part of the coefficients of the gradient 
and nonlinear terms vanishes in the long wave-length 


In this paper, we have studied nonequilibrium steady 
states of a number of driven-dissipative systems exhibit¬ 
ing a nontrivial competition between drive and dissipa¬ 
tion. Mean field theory is often used to predict the many- 
body phases and phase transitions of such systems. How¬ 
ever, a more careful field-theoretic treatment based on 
the Keldysh formalism may invalidate certain predictions 
of mean field analysis. We saw that, for example, bista¬ 
bility is an artifact of the mean field theory (model in 
Sec. HIB). Sufficiently strong dissipation may also make 
certain phases inaccessible (model in Sec. HI G in d = 2 
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dimensions), or may turn a continuous transition into a 
first-order phase transition (model in Sec. Ill A). More 
generally, the path-integral approach and even its classi¬ 
cal (saddle-point) approximation produces better results 
than mean field theory not only in equilibrium [91], but 
also away from equilibrium (model in Sec. HID). 

In all cases, an effective temperature emerges as the 
result of dissipation, and the universal behavior includ¬ 
ing the dynamics near the steady state is generically 
described by a thermodynamic universality class. The 
emergent thermal character of driven-dissipative systems 
may be expected as the quantum coherence is lost to dis¬ 
sipation. However, the phase diagram and the nature 
of phase transitions, and the precise equivalence with a 
particular thermal model is often nontrivial, and requires 
a rather careful treatment based on the Keldysh formal¬ 
ism. This paper offers such a systematic study of four 
models in great detail and, therefore, illuminates path¬ 
ways for the beyond-mean-held study of a wide range of 
other driven-dissipative systems. 

We conclude by mentioning other examples of driven- 
dissipative systems that should be amenable to a similar 
treatment. Notable models, also of experimental rele¬ 
vance, are systems with bosons, or photons, coupled to 
(pseudo-)spins on a lattice. While the two species of 
helds make the held-theory treatment more complicated, 
the photonic part is usually quadratic and can be inte¬ 
grated out at the level of the Keldysh action. The resul¬ 
tant effective model is also local due to the dissipative 
gap of the typically lossy photons. Examples of exper¬ 
imentally accessible systems of this kind include super¬ 
conducting circuits [9] , spin-boson networks [92] , strongly 
interacting Rydberg polaritons [16-19, 81], and internal 
states of ions coupled to their motion [93, 94]. We also 
remark that models closely related to that of Sec. IIIB 
arise when atoms are coupled to the electromagnetic vac¬ 
uum, and the latter is eliminated in the Born-Markov 
approximation [95], a setup that can also be accessed in 
reduced dimensions [96, 97]. However, in certain spin- 
boson-coupled systems, the Born-Markov approximation 
may not apply, but dissipation caused by external baths 
may act directly on the bosons and/or the spins. In 
this case too, the Keldysh formalism should apply. It 
would also be interesting to explore the applicability of 
the Keldysh formalism to situations involving dark states 
[24, 98] and situations involving transport, as particles 
or photons continuously enter the system at a boundary 
[16, 18, 90, 99]. 
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Appendix A: Mean field for Model III A 

A uniform ansatz for the Model in Sec. HI A yields the 
MF equations 

A = -Ay - -X, 

2 ’ 

Y = AX + J XZ -^Y, (Al) 

Z = -JXY-T{l + Z), 

where X = (a^), etc. We have defined A = 2A, J = 23 J, 
and 3 = 2d as the coordination number. The MF predicts 
a continuous transition at 

P^ - 4JA -f 4A^ = 0, (A2) 

consistent with r = 0 in the same subsection. 


Appendix B: Mean field for Model III B 


A uniform ansatz for the Model in Sec. IIIB yields the 
MF equations 

X = -AY + JYZ -^X, 

Y = AX -ClZ - J XZ -^Y, (Bl) 

Z = ilY-T{1 + Z). 

We have defined A = 2A, D = 2D, and J = 3J. Casting 
in terms of n = (1 -I- Z)/2, we find for the steady state 


4(A -f J - 2J nf + 20^ -f P^ 


n = 


(B2) 


which is similar to the mean-field equation in the 
continuum, Eq. (27). The MF equation above exhibits a 
continuous transition from a stable uniform phase to a 
bistable region with two stable uniform phases. 


Appendix C: Mean field for Model III C 

The MF for the model in Sec. HIC is derived in 
Ref. [29]. A two-site ansatz on the checkerboard sub¬ 
lattices A and B yields the mean field equations 

Aa = -2JZaYb - ^Xa, 

Ya = -2JZaXb - ^Ya, (Cl) 

= 2J{YaXb - XaYb) - r(i + Za), 

and a similar set of equations for A o H. The MF pre¬ 
dicts a continuous phase transition from a paramagnetic 
state with Y = — 1 to a staggered XY phase where spins 
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on the two sublattices are at angles 6 and —9 with re¬ 
spect to the X = y line on the Bloch sphere, see Ref. [29] 
and the figure therein. The phase transition occurs at 

J=\, (C2) 

consistent with setting r = 0 in Sec. IIIC. 


where J = 23 J and A = 2A. In the steady state, n varies 
between 0 and 1 depending on decay rates; however, the 
MF always gives X = Y = 0. 


Appendix D: Mean field for Model III D 


A uniform ansatz for the Model in Sec. HID yields the 
MF equations [n = (1 -I- (cr^))/2] 


A = - 


J (2n — 1) -I- A 


Y - 


r + r„ 


+ Kn ] X 


Y = 


J(2n-1) + A 


X- 


r + r„ 


+ Kn]Y, (Dl) 


= —Fn -I- Fp (1 — n) — . 
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